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Exact spectral truncations of the unforced, inviscid Burgers-Hopf equation are 
Hamiltonian systems with many degrees of freedom which exhibit intrinsic stochas- 
ticity and coherent scaling behavior. For this reason recent studies have em- 
ployed these systems as prototypes to test stochastic mode reduction strategies. 
In the present paper the Burgers-Hopf dynamics truncated to n Fourier modes 
is treated by a new statistical model reduction technique, and a closed system 
of evolution equations for the mean values of the in lowest modes is derived for 
m <S n. In the reduced model the ra-mode macrostates are associated with trial 
probability densities on the phase space of the n-mode microstates, and a cost 
functional is introduced to quantify the lack of fit of paths of these densities to 
the Liouville equation. The best-fit macrodynamics is obtained by minimizing 
the cost functional over paths, and the equations governing the closure are then 
derived from Hamilton- Jacobi theory. The resulting reduced equations have a 
fractional diffusion and modified nonlinear interactions, and the explicit form of 
both are determined up to a single closure parameter. The accuracy and range 
of validity of this nonequilibrium closure is assessed by comparison against di- 
rect numerical simulations of statistical ensembles, and the predicted behavior is 
found to be well represented by the reduced equations. 
© 2000 Wiley Periodicals, Inc. 
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1 Introduction 



Throughout the applied mathematical sciences one meets complex nonlinear 
dynamical systems having large finite dimension, or infinite dimension, whose typ- 
ical solutions exhibit chaotic or turbulent behavior. Often the same systems also 
reveal some organized structures that persist within the disordered fluctuations in- 
trinsic to the dynamics. One is then faced with the challenge of describing such 
systems by appropriate model equations that synthesize the coherent behavior and 
the fluctuations, and that properly account for the interactions between them. A 
common approach is first to devise some low-dimensional governing equations for 
the coherent behavior and then to add some noise, or stochastic forcing, to capture 
the effect of the fluctuations. But if the original system is already described by an 
accepted set of governing equations that are deterministic, the question arises: in 
what sense is the stochastic model consistent with the given complex dynamics? To 
address this question of principle, it is necessary to derive the proposed stochas- 
tic model from the underlying deterministic dynamics by a systematic reduction 
technique EKHIlllliaiaa. 

This general problem of model reduction for complex dynamical systems, usu- 
ally governed by nonlinear ODEs or PDEs, is conceptually identical to the generic 
problem of nonequilibrium statistical mechanics, which endeavors to connect macro- 
scopic descriptions of systems consisting of a huge number of interacting particles 
or modes to the microscopic dynamics of their constituents ll3l [T5l[T6*ll22ll32ll33ll . 
The standard approach in the field of statistical mechanics is to invoke projection 
operators that map the full phase space onto the span of some relevant, or resolved, 
dynamical variables, which define a macroscopic state. This method is used in the 
well-known Mori-Zwanzig formalism, which applies to an arbitrary set of resolved 
variables for a general Hamiltonian system ||3l[33j|7j|8l. Typically these derivations 
make use of near-equilibrium assumptions to define tractable projection operators, 
and they lead to stochastic integro-differential equations having memory kernels 
and non-Markovian noise. While the fluctuation-dissipation theorem determines a 
key relation between the autocorrelation of the noise and the memory kernel, the 
evaluation of the kernel itself is determined by ensemble averages of the certain 
derived observables propagated under a complementary projection of the micrody- 
namics. Consequently, it is normally imperative to make further approximations 
and simplifications, most notably taking a Markovian limit, to arrive at a useful 
and manageable model. As a result, systematic reductions of this kind are mostly 
confined to dynamical systems having special structures dictated by the methodol- 



ogy- 
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In recent work [31] one of the authors has proposed a new approach to statisti- 
cal model reduction, which applies to any Hamiltonian dynamics and any suitable 
choice of resolved variables, and which leads directly from the given deterministic 
dynamics to closed reduced equations without recourse to intermediate stochastic 
approximations. This method of model reduction relies on a statistical optimization 
principle. First, a parametric statistical model is adopted consisting of trial prob- 
ability densities on phase space associated with the mean values of the selected 
resolved variables. Then the reduced dynamics are determined by minimizing a 
cost functional over paths in the statistical parameter space. In this way a unique 
statistical closure is obtained once an appropriate cost functional is specified. The 
cost functional is designed to quantify the lack-of-fit of the path of trial densities 
to the underlying dynamics, and information-theoretic considerations suggest its 
form as a time-integrated, weighted, squared-norm on the residual of the trial den- 
sities with respect to the Liouville equation. A predicted, or estimated, evolution 
of the resolved variables therefore corresponds to a path in the statistical parameter 
space that best fits the underlying Hamiltonian dynamics in a sense quantified by 
the cost functional. All adjustable parameters in the closure appear as weights in 
the cost functional. These weights, which encode the influence of the unresolved 
fluctuations on the mean resolved evolution, are determined empirically. 

The closed reduced equations that result from this dynamical optimization pro- 
cedure have the desirable feature that they take the generic form of governing equa- 
tions for nonequilibrium thermodynamics ifTOll . Namely, they are systems of first- 
order differential equations in the mean macrostate, and they contain a reversible 
term, which is a generalized Hamiltonian vector field, and an irreversible term, 
which is a gradient vector field [29]. Moreover, the structured derivation of these 
equations implies that the reversible and irreversible terms in the reduced equations 
correspond to terms in the cost function that quantify the resolved and unresolved 
components of the Liouville residual, respectively. 

The purpose of the present paper is to illustrate and test this statistical reduc- 
tion method by applying it to the Truncated Burgers-Hopf (TBH) equation. The 
untruncated version of this partial differential equation has long been a prototype 
for understanding fluid shock formation among other phenomena GTTl . Its rc-mode 
Fourier-Galerkin truncation, without forcing or viscosity, and for sufficiently large 
n, has been proposed as a good testbed for model reduction strategies, since it arises 
as the truncation of a familiar, one-dimensional nonlinear PDE and yet shares the 
statistical properties of many important complex dynamical systems l26ll27l [Tll24ll. 
In particular, the TBH dynamics is chaotic and numerical simulations show it to be 
ergodic and mixing. The sharp truncation at wavenumber n also leads to backscat- 
ter of energy from high to low wavenumbers, which produces coherent behavior in 
the low modes. It has a (noncanonical) Hamiltonian structure, and its equilibrium 
distribution is given by a Gibbs measure. In addition, it has the important statisti- 
cal property that the decorrelation time of Fourier modes varies directly with their 
spatial scale ||2"6l , a property shared by important fluid dynamical systems such as 
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those relevant to geophysical applications. The expository article [27 ] describes the 
properties of the model both analytically and numerically, and the paper [ 1 ] details 
the Hamiltonian structure, the invariants and the ensemble statistical behavior. 

The abovementioned relationship between decorrelation time scale and wavenum- 
ber suggests a natural coarse-grained description in which the resolved modes are 
the lowest m modes, for m <^n. But, since every mode k for m + 1 < k < n is un- 
resolved in this reduction, the separation of time scales between the resolved and 
unresolved modes is not wide. Consequently, the TBH dynamics poses a stringent 
test of any closure methodology, and the efficacy of a closure for the TBH system 
has implications for the the design of effective reduced descriptions of many high- 
dimensional, complex dynamics systems encountered in the physical sciences. 

Our main result is an explicit m-mode reduced dynamics which has a structure 
similar to m-mode TBH dynamics itself, but with modified nonlinear interactions 
and strong, mode-dependent dissipation. Interestingly, the dissipation derived from 
the reduction procedure is not a standard diffusion, but rather a fractional diffusion. 
Furthermore, our statistical closure on m resolved modes only depends on a single 
adjustable parameter, which scales the magnitude of the dissipation. By comparing 
the solutions of the reduced equations against ensemble averages of direct numer- 
ical simulations for n = 50 and m = 5, we show that the reduced equations of 
the statistical closure reproduce the true statistical evolution across the m-mode 
resolved spectrum, and for large as well as small initial conditions. 



2 TBH dynamics and statistics 

The unforced, inviscid Burgers-Hopf equation is the PDE 

(2.1) | + |:(H =0 

for a scalar unknown u = u(x,t). This equation has been extensively investigated 
as a simple prototype for the nonlinear behavior typical of the multi-dimensional 
equations of hydrodynamics, especially the formation of shock discontinuities and 
development of turbulence. For these considerations it is natural to consider the 
Burgers-Hopf equation with a finite viscosity and with forcing, which may be de- 
terministic or stochastic. By contrast, the present work focusses exclusively on 



the free dynamics of spectral (or Galerkin) truncations the inviscid equation (2.1 ). 
These finite-dimensional dynamical system have been proposed as useful test sys- 
tems for coarse-graining strategies and model reduction methodologies, which are 
intended to be applied to much more complex dynamic systems ll26ll27l [Tl l24l . 

We study the dynamics of the projection of u onto the first n modes of the 
Fourier series for u. For 2% -periodic functions u(x,t), x € [0,2n), the truncated 
series is defined by the projection operator, 

« 1 [-271 

(P n u)(x,t) = Y z k (t)e ikx , Zk (t) = — u(x,t)e- ikx dx. 

k=-n 271 J 
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The mean z.o = Sq K u(x,t)dx is a conserved quantity under (2.1 1, and so it is set 
to throughout. For a real solution u(x,t), z-k = z\, for k = 1, . . . ,n. [z* denotes 
complex conjugate of z.] The state space for the truncated dynamics with n modes 
is therefore C", and a generic state is a point z = (zi, • . . ,z n ) G C. The govern- 
ing dynamics for u n (x,t) = (P n u)(x,t) is determined by the Galerkin truncation of 
( |2.1| ), namely, 

(2.2) ^ + ^( ? '")= - 

In terms of the Fourier coefficients, z = z(t), the Truncated Burger-Hopf (TBH) 
dynamics is accordingly the following system of n quadratically nonlinear ODEs: 

(2.3) "§7 + T L Zk t Zk 2 =0, (k= 1, •■•,«)■ 

Z k { +k 2 =k 

where k\,k,2 run over {±1, . . . , ±n}. 



The dynamics (2.3 ) exactly conserves the following two quantities: 



1 l' 2n 1 " 

(2.4) £ = — / u 2 dx= - £ N 2 , 

^■/o 2^ n 



1 /" 2?c 1 
(2.5) H = — ^ u 3 dx=- ZhZhzZk 



D i 1 +fo+i 3 =0 



The quadratic invariant, defines the total energy of the truncated system (2.3). 



Curiously, it is the cubic invariant, H, not the energy, E, that appears in the Hamil- 



tonian form of the TBH dynamics (2.3 1; that is, 



-77 = Y.JkW^ri; 1 f° r hk! = -ikSkk' , 
dt f dz k 

where 8kk> = 1 if k = k', and = if k / k'. These equations have the form of a 
noncanonical Hamiltonian system, or Poisson system, with cosymplectic matrix 
/ and Hamiltonian H, having n complex degrees of freedom. Nonetheless, the 
invariant H plays a much less prominent role in conditioning the dynamics than 
does the energy E, as has been documented by previous researchers ETl fTll. 

Except for n = 1,2 the TBH dynamics is chaotic, and for n > 20 numerical 
simulations indicate that it is ergodic and mixing. We therefore adopt a statisti- 
cal mechanical perspective that focuses on the evolution of ensembles of solutions 
of ( 2.3[ ), rather than individual solution trajectories. The Hamiltonian structure 
of ( 2.3[ ) implies the Liouville property, namely, the invariance of phase volume, 
dz = dx\dy\ . . .dx n dy„, (zk = Xk + %) under the phase flow. Consequently, the 
propagation of probability by the TBH dynamics is represented in terms of a proba- 
bility density p(z,t) with respect to dz that is transported according to the Liouville 
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equation, 

(2 6) o = dp i y ± k dp i y k dp = dp i y z dp 



In principle, p(z,t) is completely determined by (2.6), given an initial density 
p(z,0), and hence the expectation of any observable F on C" at time f > is deter- 
mined by 



(F(t)\p) = [ F(z)p{z,t)dz. 

JC" 



[Throughout the paper (F\p) denotes expectation.] In practice, however, the ex- 
act density p(z,t) evolves under the nonlinear, chaotic dynamics in an extremely 
complicated way, and numerical evaluation of the expectation of any observable 



with respect to p(z,t) requires computing an ensemble of trajectories (2.3 1 starting 
from a sufficiently large sample of the initial density. This intrinsic feature of the 
dynamics, which is shared by many turbulent dynamical systems, is the fundamen- 
tal reason for invoking model reduction procedures, which are devised to furnish 
sufficiently accurate approximations to those expectations without expensive sam- 
pling of the detailed dynamics. 

Our reduction technique is based on using a convenient family of trial proba- 
bility densities on the phase space to approximate the exact density evolving under 
( |2.6| ). This family of nonequilibrium densities is constructed from a fixed reference 
equilibrium density, which we take to be the canonical ensemble, 



exp(-0£(z)) i/3 _ 



<2J) p " iz) = i,M-mz))dz = U" 



e 



* * 



with specified inverse temperature j3 > 0. We choose the first m components, 
A(z) = (zi, • • • ,z m ), of the microstate z £ C" to be the resolved variables, or re- 
solved modes, in the model. Experience with computed solutions of the TBH 
equation, as well as a heuristic scaling argument, informs this choice, because it is 
observed that the decorrelation time scales for the modes scale with \/k E81I231 . 
Accordingly, the lowest m modes can be considered slow variables in the n-mode 
TBH dynamics, and the resolved vector A constitutes a natural coarse-graining 
of the microstate z £ C. We call the mean values A the macrostate and denote 
it by a = (ai, . . . ,a m ). This coarse-grained description choice does not, however, 
ensure a wide separation of time scales between the resolved and unresolved vari- 
ables, and consequently the TBH dynamics presents a rather challenging test for 
the model reduction technique. 

The general properties of the non-equilibrium probability distributions of this 
system were examined numerically in lTT8l where it was noted that they are of- 
ten quasi-normal. Motivated by this empirical observation, our statistical model 
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employs the quasi-equilibrium (or quasi-canonical) densities 
(2.8) p(z;X) = exp 




This exponential family of densities ( |2.8[ ) form a parametric statistical model, the 
parameter vector being X = [X\, . . . , X m ) S C m (51 [191. The convention, = X£, 
fork = 1, . . . ,n, and Ao = 0, applies to the parameters. Since pp is a Gaussian 
density, each density p(X) is also Gaussian; with respect to p{X) the means and 
variances of the resolved variables are 

(2.9) a k = (zk\p(X)) = pX k , {{zk-a k ){zie -a v )*\p(X)) = -^8 kV , 



for kjk 1 = ±1, . . ., ±m. A highly desirable feature of the statistical model (2.8 ) is 
that all the modes Zk are independent and the correspondence between the mean 
resolved variables a k and the statistical parameters X^ is linear and decoupled. 

Other possible choices of statistical model are conceivable. For instance, the 
invariance of both E and H suggests using the more general equilibrium density 

Pp,a, v {z) = Q((3,a,Ti)- l exp(-pE(z)-r,E(z) 2 -aH(z)) 

with tj > 0; <2(/3 , a, T]) is the normalizing factor. This density is a so-called Gauss- 
ian ensemble with respect to E, as it includes the quadratic term E 2 in the expo- 
nential, which is necessary to make the density normalizable given the presence of 
the cubic invariant H. But, using pp.a.r) as the reference density for the statistical 
model p(X) leads to non-Gaussian trial densities, and hence an intractable sub- 
sequent analysis. Fortunately, simulations of ensembles of solutions of the TBH 
dynamics have shown that the distribution of the modes Zk is very nearly Gauss- 
ian, and that H has only a slight effect on the statistical behavior, provided that the 
mean value (H) is close enough to the relaxed value (H\pp) =0 ir2"8ll2"3l lTI. Ac- 
cordingly, our statistical model and closure procedure are built from the reference 
density ( |2.7[ ) that does not include H. In this formulation we expect to approximate 
a statistical evolution in which the mean value of H is near its relaxed value, but 
not an evolution from an initial statistical state with a large mean value of H. 

Neither (E\p(X)) nor (H\p(X)) is exactly conserved along an arbitrary ad- 
missible path X = X(t) in the statistical parameter space. It is possible to al- 
ter the formulation of the statistical model to ensure that both of these invariants 
are respected, but with the inevitable consequence that the trial densities are non- 
Gaussian [31 ]. Our main justification for allowing some variation in these dynam- 
ical invariants is, therefore, expediency, combined with the fact that for n S> m the 
relative variations in the means of E and H are small. 

3 Formulating the optimization principle 

Our m-mode closure of the n-mode TBH dynamics is derived from an opti- 
mization principle over paths, X(t), in the parameter space, C"\ of the statistical 
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model (2.8 ). A predicted, or estimated, evolution of the macrostate a = a(t), with 
ak(t) = {zk \p(X(t)) ) for k = 1, . . . ,m, corresponds to a path X(t) that is optimally 
compatible with the underlying TBH dynamics, in the sense that X{t) minimizes 
a certain lack-of-fit cost functional over all admissible paths. This lack-of-fit is a 
metric on the Liouville residual, a fundamental statistic in our approach that we 
define to be 

\ k=—n " / k=—m 



Here dzk/dt is given by the governing TBH dynamics (2.3 1, and X k = dX k /dt. 



If p(- ; X{t)) is replaced by an exact solution p(- ,?) of (2.6) in (3.1 1, then R = 
identically on the phase space C". In OTI it is shown that the statistic R = R(z; X , X ) 
may be interpreted as the local rate of information loss at the sample point z due to 



reduction via the statistical model (2.8 1. The significance of R is also revealed by 
the family of identities 

j t (F\p(X(t)))-(F\p(X(t))) = (FR\p(X(t))), 

which hold for every observable F on C", with dF/dt = 0; here, F = {F,H} de- 
notes the Poisson bracket of F with H, so that dF(z(t)) jdt = F(z(t)) on solutions 



z. = z.{t) of (2.3 1. Since the mean of R with respect to p(X) is zero, (R\p) = 0, 
the covariance between an observable F and the Liouville residual R quantifies the 
deficiency of the path of trial densities p(X(t)) to propagate the expectation of F. 
This family of identities over test observables furnishes a natural linear structure 
for defining the lack-of-fit of any admissible path X (t) to the underlying dynamics, 
and motivates the use of a weighted L 2 (C n , p (X ) ) squared norm as the lack-of-fit 
cost function. 

The weights in the lack-of-fit cost function relate to the resolved and unresolved 
subspaces of L 2 (C n ,p(X)). The mean-centered resolved variables 

(3.2) U k = Zk-a k = 3 ^ r logp(-;A) (k=l,...,m), 

oX k 

span the resolved subspace, and, by ( |2.9| >, they are uncorrected Gaussian variables, 
(U k U^,) =fi- x 8 kk! , for all XeC". The orthogonal projection of L 2 (C",p) onto the 
resolved subspace is 

m 

(3.3) P A F = P £ (FU* k \p)U k , 

k=—m 

and the complementary projection onto the unresolved subspace (the orthogonal 
complement of the resolved subspace) is Qa = I — Pa- The projections of the 
Liouville residual R = R(X,X) onto its resolved and unresolved components are 
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(3.4) 



p a r = p £ 



k=—m 



ik 

dk+ 2 



a ky a kl 



ki +k 2 =k 



(3.5) Q A R 



ka *k L U h u k 2 = l ^ £ k 3ak 3 U kl U k2 . 



k=—m k\+k 2 =k 



ki+k 2 +k 3 =0 



Here, d k = da k /dt = X k /f5. Since a k = for \k\ > ra, the wavenumber indices k\,k2 



in (3.4 1 run over {±1, . . . , ±m}, unlike the nonlinear term in the TBH dynamics, 
in which all modes in {±1, ... , ±«} interact with the resolved modes. By contrast, 



£i,&2 in (3.5 1 run over {±1, . . . , ±n}, and the definition of U k is extended for \k\ > 



m to be simply U k = z.k- 



The identities ( 3.4 ) and ( 3.5 1 make use of the following identities for the Gauss- 
ian scores variables: 



(3-6) 



(U kl U k2 U k3 \p) 



(U h U k2 U^\p) 



0. 
1 



The calculations (3.4 1 and (3.5) are summarized as follows. In light of (3.1 ) 



dzk 



p a r= £ X;u k +W £ (-£u})u k 



k=—m 



k'=—m 



Substituting (2.3 1 and using P(z kl Zk 2 U k , ) = a kl 8 k2k > + a kl 8 k{k i , which is implied by 



(3.6 1, results in 



P A R= £ X* k U k + X* k ik £ a^uUtf 

k=—m k'=—m 

m | m | 

= £ \K-i £ m*u\u k 

k=—m I k'=—m I 



£ + ' £ hX kl a k2 > U k . 

k=—m k\+k 2 =k ) 
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Putting Xk = ficik and using the symmetry between ki and &2 then produces the 
desired identity (3.4). In a similar manner, 

Q A R = R- P A R 

dzk 



22 K 

k=—m 



dt 



ik 22 

k'=-r. 



O-k-k' 



22 

k=—m 



kX 



■k 22 Zk l zk i 

k[+k 2 =k 



ak,Uk 2 —ak 2 Uk, 



22 *K 



* I 

k=—m k[+k 2 = 



Uk { Uk 2 —ak t ak 2 



using the symmetry between k\ and &2 in the third equality. The constant term in 
the last expression vanishes, because 



k=- 



22 a h a h 

k\+k 2 =k 



~P 22 k^a^a^a^ 

ki+k 2 +k 3 =0 



0, 



and hence the desired identities (3.5 ) follow. 



The cost function for our optimization principle is declared to be 



(3.7) J?(X,X) = -([P A R(X,X)] Z \m 



+ l 2 ([WQ A R(X,X)] 2 \p(Xy, 



in which W is a linear operator on L 2 (C",p) satisfying WP A = P A W. We refer to 
W as the weight operator, because it weights the contributions of the unresolved 
component, Q A R, of the Liouville residual R relative to the (unit weighted) resolved 
component, P A R. The requirement that W commutes with the projection P A implies 
that W also commutes with Q A and hence that W takes the unresolved subspace into 
itself. The inclusion of W in the cost function (3.7 1 gives our best-fit closure the 
character of a weighted least-squares approximation over paths X(t). 

With this cost function in hand, we are now able to formulate the optimization 
principle that defines our statistical closure. In this section we present the stationary 
version of the principle, which is simpler to describe and motivate. A nonstationary 
version, which is needed when comparing the predictions of the closure with direct 
numerical simulations, is given in Section 6. 

The best-fit closure scheme is based on the dynamical minimization problem 



(3.8) 



v(A°) 



mm 

A((o)=A° 



Sf(X,k)dt, 



in which the admissible paths X(t), to < t < +°°, in the configuration space of 
the statistical model are constrained to start at X° G C at time to. In optimization 
and c ontrol theory, v(X ) is called the value function for the minimization problem 
(3.8 1 [i4j[TTl[T2l. Since j£f (X, X) is independent of t, and the integration extends to 
infinity in time, v(X) is time-independent, and the initial time to may be shifted to 
Oinjia. 
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By analogy to analytical mechanics, one may regard (3.8 1 as a principle of le ast 
action for the "Lagrangian" Jzf(A,A) and interpret the first member in (4.6) as 
its "kinetic" term and the second member as its "potential" term |f2j[T3j|20l. The 
kinetic term is a positive-definite quadratic form in the generalized velocities X , and 
it is entirely determined by the structure of the resolved variables and trial densities. 
By contrast, the potential term embodies the influence of the unresolved variables 
on the resolved variables and involves the weight operator W, which contains all 
the adjustable closure parameters. 

Each extremal path X(t), to < t < +°°, for (3.8 1 corresponds to an evolving 
trial density p( • ; % (/)) that is best-fit to the Liouville equation in the sense that the 
time-integral of the lack-of-fit cost function is minimized. Finiteness of the value 
function implies that X(t) — > as t — > +°°, meaning that the best-fit path connects 
the given initial state A to equilibrium X eq = 0. Hence, these extremal paths model 
the relaxation of the mean resolved vector d(t) = (A \p(X(t))), from a given initial 
stated = (A|p(A )). 

It is important to emphasize that the weight operator W is an independent ingre- 
dient in our closure strategy and that W includes all the empirical parameters in the 
best-fit closure. It is unavoidable that some empirical ingredient should enter into 
the definition of such a closure, because the defining optimization principle does 
not refer to the autocorrelations of the unresolved modes under the full dynamics, 
or the projected dynamics orthogonal to the resolved subspace. As is well-known 
in nonequilibrium statistical mechanics, the transport properties of the complex 
system, as well as the memory kernels in its projected dynamics, are expressible 
in terms of such autocorrelations 121 |32l [33]]. In this light our closure strategy 
of adopting a time-integrated, weighted, least-squares approximation is a practi- 
cal expedient to avoid the expensive computation of such autocorrelations, which 
require the propagation of ensembles of solutions of the fully resolved dynamics. 
Our approach instead quantifies the influence of the unresolved modes on the re- 
solved modes by an appropriately chosen weight operator W which depends on 
some adjustable parameters that must be tuned empirically. 

In this general formulation of the optimization principle, the weight operator W 
is not restricted to any particular form. Moreover, the closed reduced equations sat- 
isfied by its extremals have a generic thermodynamic format and associated prop- 
erties for any admissible choice of W OTTl . In any concrete implementation of the 
best-fit approach, however, it is essential to discern a convenient form for W from 
the structure of the underlying dynamics, the resolved variables and the unresolved 
Liouville residual. It is also desirable to choose a weight operator that contains as 
few adjustable parameters as possible, given the required level of fidelity of the clo- 
sure approximation to the true statistical dynamics. In the Burgers-Hopf dynamics 
the quadratic nonlinearity, typical of hydrodynamics, affords a particularly simple 
and efficacious choice of W, as is described in the next section. 
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4 Deriving the closed reduced dynamics 

First, let us consider the implications of the trivial choice for the weight op- 
erator in ( |3.7[ ), namely, W = 0. In that situation the solution of the minimization 
problem ( |3.8[ ) is complete determined by setting = at each instant of time. 
The resolved vector, a = (a\(t), . .. ,a m (t)) = ((zi\p), ■ ■ . , (z m \p)), then satisfies 

da k ik 
— - + — 
dt 2 



+ T L a h a h =° (k=l,...,m) 



meaning that the reduced dynamics coincides with the spectral truncation to m 
modes of the Burgers-Hopf equation. Not only are the m-mode truncations of E 
and H invariant under this dynamics, but the entropy, 

o m 

(4-1) s = -(logplp) = ~ E N 2 , 

m 

is also invariant; indeed, the entropy production is 
ds ^ ^da k ip 

al k=-m al Z ki+k 2 +h=0 

This naive closure is therefore adiabatic, in that it suppresses interactions between 
the resolved modes and the unresolved modes, those interactions being the source 
of entropy production and dissipation in a proper reduced dynamics. 

Now, we seek an effective choice of W that appropriately quantifies the cost of 
the unresolved component of the Liouville residual as well as the resolved com- 



ponent. An examination of (3.5 1 reveals that unresolved residual, QaR, is a linear 
combination of the products U kx U k2 over &i,&2 = ±l,...,±n. This structure is 
a consequence of the quadratic nonlinearity of the TBH dynamics. Since these 
products form an orthogonal basis for the unresolved subspace, it is natural and 
convenient to construct the desired weight operator in terms of them. Accordingly, 
we take W to be diagonal with respect to this basis and we set 

(4.2) W(U kl U k2 ) = s klM (U kl U kl ), 

for some sequence of constants e kl jt 2 > 0, symmetric with respect to interchanging 
&i with ki. The "potential" term in the cost function (|3.7|) is then 



1 1 m 

-([WQ A R(X,X)] 2 \m) = ^ J Lk 2 \a k \ 2 - £ el M (U h U k2 U^ 2 ) 

k=—m k[+k 2 =k 
i m 

(4-3) Jkk 2 \a k \\ 

Z P k=-m 

with 

(4-4) y k =\ £ t\ M (k=l,...,m). 

£ kl +k 2 =k 
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The sums in (4.3 1 and (4.4 1 extend over 1 < \ki\, I&2I < n. The calculation in (4.3 1 



uses the fourth moments in (3.6 1. We thus arrive at the explicit expression for the 
cost function, 



1 m 1 

(4.5) j2f(A,A) = 2 £ « 



Z P k,+k 2 =k 



+ tt*ia 1 2 



which follows from (|3.7|) by (|3.4|) and (4.3 1, together with the simple relations 



given in ( |2.9[ ). 

We notice that the weights £k u k 2 defining W enter the cost function Jz? only 
though the combinations 71 , . . . , J m in (|4.4|). While it is possible to retain these 



m adjustable parameters in the closure, their form as sums over the weights £k u k 2 



suggests a further simplification. Namely, for 1 < k < m <C n, the % in ( 4.4 ) are 
approximately independent of k, whenever the weight sequence ^ is chosen to 
be bounded and slowly varying over the range 1 < \k\\, I&2I < n. On this basis it 
is reasonable to set 7. = 7 > 0, for k = 1 , . . . , m, and to adopt a cost function with 
a single adjustable parameter, 7. We confine ourselves to this simplified formula- 
tion throughout the present paper. This choice of cost function is also justified a 
posteriori by our subsequent analysis, which shows that the scaling of the modal 
dissipation rates with \k\ implied by a constant 7 agrees with the direct numerical 
simulations presented in Sections 7 and 8 as well as the formal dimensional ar- 
gument summarized in Section 5. Indeed, we find that the TBH dynamics is an 
example of a complex system in which the effect of unresolved modal fluctuations 
on the resolved modes is efficiently and accurately modeled by a lack-of-fit cost 
function having a single, empirically-tuned, closure parameter. 

The equations satisfied by the extremal paths for ( |3.8| > are found by the tech- 
niques of the calculus of variations ||4l [T21 [HI El [30l . In particular, Hamilton- 
Jacobi theory furnishes the the closed reduced equations for the statistical model 



in terms of the value function v(A). To this end we write the cost-function (4.5 1 in 
the more compact form, 



(4.6) 



m 1 
k=l P 



h-pfk(X) 



+ rk 111 2 



introducing the shorthand notation 



(4.7) 



fk(X) 



ik 



ki+k 2 =k 



XkyXk 2 . 
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and recalling that A_£ = X£. We form the Hamiltonian J4?(X,jJ.) conjugate to the 
Lagrangian ££ (A, A) by taking the Legendre transform of Jzf: 

<9^f 1 • 

(4-8) ^i k = = X k -f k (X), 

m 

J$?(X,ii) = ^^ k X k + X^ k -^(X,X) 
k=\ 

m y 

= Y,m 2 +^Mfk{x)+f k {xr P i k \-j i k 2 \x k \ 2 . 



k=\ 



[That these complex expressions are appropriate complexifications of the real Le- 
gendre transform is shown in the Appendix.] 



The value function, v(X), in (3.8 1, which is analogous to an action integral or 



Hamilton principal function, satisfies the stationary Hamilton- Jacobi equation: 
(4-9) ^(A,-^)=0, 

According to Hamilton- Jacobi theory, the conjugate variable /I = jX(t) along an 
extremal path A = A (?) is given by the relation 

(4.10) f = 

This key relation closes the reduced dynamics along the extremal path, since to- 



gether with (4.8 1 it yields the system of first-order ODEs 



(4-11) d ^r =fk(X) -^ 1 ^ (k = l,...,m) 

recalling that a k = X k /f5. When expressed entirely in terms of the mean resolved 
variables, the closed reduced equations are 

(4.12) ¥ +2 t|+ I%^ = -^(H. 

The left-hand side of this system of equations is recognizable as the m-mode 
TBH dynamics. The adiabatic closure mentioned earlier in this section is obtained 
by setting v = identically in ( 4.12| ). The presence of the gradient vector field of 



v on the right hand side provides the dissipation, or irreversibility, of the closed 
reduced dynamics. Indeed, the relation ( |4. 10 > has the interpretation that the conju- 



gate vector /I to the statistical state vector A represents the irreversible part of the 
flux of the mean resolved vector a. 

The best-fit closure ( 4.11 14. 12 > is completed once the value function (3.8) is 



determined explicitly, at least to some suitable approximation. This calculation is 
the content of the next section. 
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5 Approximating the value function 

Before presenting the nonlinear analysis needed to approximate the value func- 
tion, we exhibit the first-order, linear approximation in the near-equilibrium regime, 
which can be obtained by elementary methods. We consider the relaxation of a 
small initial disturbance from equilibrium, and so we retain only the quadratic 
terms in the cost function 04.6b, which becomes 



^(A,A) = £i|A,| 2 + ^|4| 2 . 
k=i P P 



Each extremal, X(t), solves the Euler-Lagrange equations for (3.8 1 for this inte- 
grand, 

Since the admissible paths connect the initial state to equilibrium as t — > +oo, the 
relevant extremals are the exponentially decaying solutions, 

X k (t) = 4(0) exp . , ; g 

The near-equilibrium closed reduced equations satisfied by these extremals are 

(5.1) § ^ 

at 

and the value function evaluated on these extremals is 

rv m 





The elementary result (5. 1 ) is consistent with ( 4.11[ >, since / = in this approx 



imation. Moreover, the same result is anticipated by a heuristic scaling argument in 
|[26ll27l . Phenomenological reasoning about the transfer of energy in the turbulent 
TBH system suggests forming an eddy turnover time, T k , for the k-th mode from 
the wavenumber, k, and the energy per mode in statistical equilibrium, 1//3. The 
only dimensionally consistent time scale is 7* ~ \/P/k. It is well confirmed by 
direct numerical simulations of the TBH dynamics in Il26ll27ll that the equilibrium 
autocorrelations of the modes, z.k, exhibit the ^-dependent decay rate proportional 
to TT 1 across the spectrum, thus validating this scaling argument. To relate our ba- 



sic prediction (5.1 1 to this scaling theory and its validating numerical simulations 
we rely on linear response theory, which holds for relaxation near equilibrium and 
connects the nonequilibrium mean value of the resolved variables A* to their equi- 
librium autocorrelations ||6l [15] |33]]. Precisely, the ensemble mean of A k (t), the 
resolved variable propagated under the phase flow for time t, with respect to an 
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initial quasi-equilibrium density p(A°) for small is calculated up to linear- 
response approximations to be 

m m 

(A k (t)\p(X°))^(A k (t)[l + £ {X»yA k ,]) eq = £ (A k m,) eq X°. 



k'=—m 



k'=—m 



Thus the time scale for decay of (A k (t) |p(A )) is dictated by the time scalings 
of the decaying autocorrelations (A k {t)A* k ,) eq . In light of the robust scaling be- 



havior of the autocorrelations, the decay rates in (5.1 ) are necessarily proportional 
to \k\/y/p. This correspondence also validates the choice of the weight operator 
in the cost function and the further simplification to a single closure parameter y, 
independent of k. 

We now turn to the nonlinear problem of determining an approximation to 
v{X) that is valid for larger amplitudes We seek an solution that is accu- 



rate up to cubic terms in A, so that the resulting gradient term in (4.12 1 contributes 
quadratically-accurate nonlinear terms. Accordingly, we submit a third-order Tay- 



lor expansion for v to the Hamilton- Jacobi equation (4.9 1 : 



(5.2) 



v(A) = £m*|A* 



k=i 



+ £ N k{klkj X k{ X kl X kl +. 



The constant and first-order terms in (5.2 1 vanish, since v(0) = = dv/dX{0). The 
quadratic terms in the expansion ( |5.2| already have the special form anticipated by 
the foregoing linear analysis, with real and positive coefficients M k . The indices 
fcl ,&2) £3 for the cubic terms run over { ± 1 , . . . , ±m} , and the complex coefficients 
Nkfah wt assumed to be symmetric under permutation of &i,&2,^3 and to satisfy 
N-ki.-k2.-k3 = A^ ]jt j , since v(X) is a real- valued function. These requirements 
supplemented the convention that X- k = Xj* and \L- k = \i k . The reader is referred 
to the Appendix for a summary of the appropriate form of Taylor's theorem for 
functions of several complex variables. 



The Hamiltonian in (4.9 1, when written in full using (4.7 1, is 



(5.3) je{X,n) = £fc| 2 -^ 2 |4| 2 + -L £ hX kl X k2 Vh- 

P Z P ki+k 2 +k3=0 



k=l 



The leading term in (4.9 1 is evaluated as follows: 



(£4) 



k=A d K 



™ dv dv 
dX k dX_ k 

m m ( d d 

V M k \X k \ 2 + V N hhh Y M k X k ^r- +X- k ^r— ) X h X kl X h 



k=l 
111 



^M^Xkl 2 + £ N kihh {M h +M kl +M h )X h X kl X h . 
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Substituting dv/dX^ = M k X k + 0(\X\ 2 ) into (4.9 1, and using (5.4 1, yields the Hamilton- 
Jacobi equation up to third order: 



(5.5) Y<P M k\h\ 2 --L£\h\ 2 + P £ N klk2k3 (M kl +M k2 +M k3 )X kl X k2 X, 



'&3 



£ hM ki X k{ X k2 X k2 = 0. 



Equating the quadratic terms in \5.5\ produces the anticipated coefficients, 

| (k = l,...,m). 



(5.6) 



M, 




Equating the cubic terms in ( |5.5[ ), and symmetrizing the second cubic term, pro- 
duces 

j3 £ N klk2k3 (M kl +M k2 +M k3 )X kl X k2 X k3 

h,k 2 ,h 



6j3 



£ (k\M h + k 2 M k2 + k 3 M k3 )X kl X kl X kl = . 



ki+k 2 +k 3 =Q 

It follows that 

N kl k 2 k 3 = , unless ki+k 2 + k 3 =0, 
and for index triples with k\ + &2 + k 3 = 0, 

i kiM kl +k 2 M k2 +k 3 M k3 



(5.7) 



N hk 2 h 



6/3 2 M^+M^+M^ 
6F |*i| + |* 2 | + |* 3 | 



in light of ( 5.6 1. We note that N klklki is symmetric under permutation of its indices, 
and that N- kl ,-fc 2 ,-£ 3 = ^ kl kl k3 - Thus, the coefficients of the approximate solution 
( |5.2[ ) are determined. 

This explicit expression for the value function v(X) furnishes the gradient term 
juatio 

M k X k 



in the closed reduced equation (4.12); namely, for k = 1, . . . ,m, 

dv ... d 



i) L N klklki X kl X k2 X k3 

0/l k fciH-fc 2 4-fc3=0 

-M k X k -3 £ N kltk2 - k X kl X k2 

kl +k 2 =k 

... i T-i k\\ki\+k 2 \k2\-k 2 



using the symmetry of N k{klk3 in the second equality. 
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The closed reduced equation, accurate to second order in amplitude, for the 
resolved vector a = ( (zi | p) , . . . , (z m I P) ) is thus found to be 

dci iJc 

dt 1 h+k 2 =k 




in which 



CO{k\,k 2 



k Y \k l \+k 2 \k 2 \-{k i +k 2 Y 



(*i+* 2 )(|*i| + \H+ k \ +h) ' 
As in the TBH dynamics itself, the indices k\ and k 2 in the convolution-like sum run 



over {±1, . . . , ±m}. The effective dynamics (5.8 1 for the m resolved modes is our 
main result for the stationary version of the best-fit closure scheme. The dissipative 
terms coincide exactly with those obtained above in the near-equilibrium, linear- 
response regime. The quadratic interaction terms in ( |5.8| ) have the same form as the 



corresponding terms in the TBH dynamics (2.3 1, but they are modified by addition 
of the factors co(ki,k2). Thus, the best-fit closure procedure produces governing 
equations for the mean resolved modes which capture both the linear dissipation 
and the modified nonlinear interactions of the resolved modes which result from 
couplings between the resolved and unresolved modes. 



The dissipative structure of ( 5.8 1 is novel in the sense that it does not correspond 
to a standard diffusion, for which the decay of mode k scales with k 2 . Rather 
the dissipation for the closure of the TBH dynamics is a fractional diffusion, for 
which the decay of mode k scales with \k\. As mentioned above, this prediction of 
our closure theory is robustly observed in numerical simulations, including those 
presented in Sections 7 and 8. 

The effects of the modification of the TBH mode interactions are more subtle. 
The factors co(ki,k2.) can be either positive or negative depending on the signs of 
ki and k 2 . Specifically, 

1 —k\k 2 
— < (0(k u k 2 ) = < 0, when ky > 0, k 2 > 0, 

4 [k\+k 2 ) z 



< co(ki,k 2 ) = — <1 when ki>0,k 2 <0, 
k\ 

and symmetrically when k\ <0,k 2 > 0. Accordingly, the low-to-high mode trans- 
fers for the closed reduced equation are weaker than for the TBH equation itself, 
while the high-to-low mode transfers are stronger. Also, there is an asymmetry be- 
tween the weakening and the strengthening of the mode transfers: The downscale 
(low-to-high mode) transfers are weakened by factors of at most 25%, while the 
upscale (high-to-low mode) transfers are strengthened by factors that can approach 
100%. 
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6 Nonstationary formulation of the closure 

The stationary formulation of the best-fit closure scheme developed in the pre- 
ceding section produces a time-independent value function, v(A), which deter- 
mines the irrev ersible component of the flux of the estimated mean resolved vector 



a according to (4.1 1 1. The entropy production along an extremal X(t) for the sta- 



tionary optimization principle ( |3.8| > is given by 

ds 2k * dv ,s > 
dt = L h dX k {X) - 

k=—m K 

It follows that throughout the subdomain over which the value function is convex, 
which includes a neighborhood of equilibrium, 

f >v(X)>0, 

since, by convexity, = v(0) > v(A) — Y, m m Xkdv/dXk. Conceptually, this inequal- 
ity is entirely satisfactory: the entropy production is bounded below by the value 
function, which quantifies the optimal information loss in the closure and depends 
only on the statistical state A. But, when a quasi-equilibrium density p(A°) is 
specified as the initial state at t = 0, the entropy production of the path of quasi- 
equilibrium states emanating from p(A ) and following the exact mean values a(t) 
tends to zero as t — > 0+. In this situation the entropy production, and the irre- 
versible component of the flux, vanish initially and there is a "plateau" time over 
which they grow to their stationary values OTTl . Since our numerical experiments 
in Sections 7 and 8 make use of such initial conditions, it is necessary to modify 
the formulation of the defining optimization principle to include this plateau effect. 
We refer to this modification as the nonstationary formulation, in that the value 
function v = v{X,t) becomes time dependent, and v(A,0) = identically. 
The nonstationary optimization principle has the value function 

(6.1) v(AVi) = min P Sf(X,-X)dt, 

A(f 1 )=A 1 Jo 

in which the admissible paths X(t), < t < t\ , in the configuration space of the 
statistical model are constrained to terminate at A 1 G C" at time t\ > 0, while X (0) is 



unconstrained. The integrand in ( 6. 1 1 is modified to account for time-reversal. The 
admissible paths may be viewed as evolving in reversed time, % = t\ — t, starting 
from X 1 at z = 0. The value function at t = t\ then quantifies the optimal lack-of-fit 
of time -reversed paths connecting the current state X 1 to some (unspecified) initial 
state X (0) corresponding to some (quasi-equilibrium) trial density at t = 0. A fuller 
justification is given in OTI . 



In terms of this value function ( 6. 1 1, closure is achieved at each time t by setting 
(now with (A,?) replacing (A \t\) ) 

(6.2) £(,) = __i!_(A,/). 
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This relation is the nonstationary analogue to (|4.10|), and it determines the reduced 



dynamics in the same way that (4. 1 1 1 follows from (4. 10 1. The resulting equations 
for the mean resolved variables are the nonstationary version of ( 4.12| ); namely, 



(6.3) 



ddh ik 

— - + — 
dt 2 



k 1 +k 2 =k 



dv 



The value function ( 6. 1 1 is the unique solution of the initial value problem 

dv 



(6.4) 



dt 



+ je\x 



dv 

d~x* 



0, for t >0, with v(A,0) =0, 



-jX ) is the Le- 



which is a time-reversed Hamilton- Jacobi equation; that is, J#°(A, 
gendre transform of Jzf (A, — A). 

Since J$?(X,n) is positive-definite in fx, the solution v(X,t) of (6.4) exists for 
all time t > 0, and as t — > +°°, v(X,t) — > v(A), the stationary value function. Thus, 
the nonstationary best-fit closure is a natural generalization of the stationary closure 
that straightforwardly includes an intrinsic plateau effect. 

As in Section 5, the closure that holds in the near equilibrium regime is calcu- 
lable by elementary methods. In the defining nonstationary optimization principle 
(6.1 1 the extremal path X(t) satisfies dX/dt(0) = 0, because A(0) is free; solving 
the Euler-Lagrange equations with this initial condition gives 



cosh , 



X k (t)=X k (h)- 



l\k\t 



cosh . 



l\k\h 



an extremal that decays exponentially in reversed time. The near-equilibrium, non- 
stationary, closed reduced equations are accordingly 



(6.5) 



dd k 
~dt 



\k\ tanh 




These equations clearly approach the corresponding stationary equations ( 5. 1 1 asymp- 
totically as t — > +oo. The value function evaluated on this family of extremals is 
calculated to be 



v(X,t) 



« 3 £*i**i 

p k=\ 



tanh 




Proceeding to the quadratically nonlinear corrections for the nonstationary closed 



reduced equations, we consider a Taylor expansion (5.2 1 for v(X,t), in which now 
the coefficients are time-dependent, M k (t), A^ Ij/ t 2 ,/t 3 (f)> andMfc(O) =0, N^^-ki (0) = 
0. The analysis of Section 5 then follows along the same lines as in the stationary 
version, except that now M k and N ku k 2 ,h satisfy ordinary differential equations 
rather than algebraic equations. Specifically, the quadratic and cubic terms in the 



expansion of the time-reversed Hamilton- Jacobi equation (6.4) yield the ODEs 

dM k 



(6.6) 



dt 



^k 2 

j3 2 : 
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(6.7) 



dN k . 



k^k 



3 -+P(M kl +M k2 +M k3 )N kl 



The solution of the Riccati differential equation (|6.6|) with M k (0) = is 



-{k\M h +k 2 M k2 +k 3 M k3 



M k {t) 



\k\ tanh 



-\k\t 



this result produces the linear terms in (6.3 1, which are anticipated in (6.5). In- 
clusion of the quadratic terms in (6.3 1 associated with the coefficient functions 
N klk2k3 (t) produces the nonstationary version of the closed reduced equations (4.12). 
Namely, for k = 1 



,.,m, 



da k ik 



£ [\+Q.(k u k 2 ,t)]a kl a kl 



\k\ tanh 



ki+k%=k 

where the functions £l{]c\,ki,i) solve the ODEs 
dQ.(k\,k2,t 



\k\t d k , 



(6.9)- 



dt 



+ p [M h (t) +M k2 {t) +M k (t)\a{k u k 2 ,t) 
= £ [hM kl {t)+k 2 M k2 {t)-kM k {t)} 



[k\ +k 2 



along with homog eneo us initial conditions Q.(k\,k 2 ,0) = 0. These ODEs follow 
immediately from (6.7 1 by setting Q.(ki,k 2 ,t) = (6fi 2 /ik)N kl}k2 .- k (t) for k\ +k 2 = 
k>0. 

Even though the equations for the factors Q.(k\,k 2: t) are inhomogeneous lin- 
ear ODEs, their coefficients and source terms are time-dependent, and hence they 
are not solvable analytically. A simple approximation, which may be helpful 
conceptually, is to replace the functions M k (t) in \6.9) by their saturated values, 
lim f _ >+0 oM ( t(f). Then the time-dependent factors £l(k\,k 2 ,t) are simply related to 
the time-independent factors 0(^1,^2) in (5.8 1; namely, 



£l(ki,k 2 ,t) & co(ki,k 2 )-< 1— exp 



|*l| + |*2|+*l+*i]f 



It is transparent from these approximate expressions that the modifications of the 
nonlinear interactions in the nonstationary formulation saturate to the stationary 
formulation with saturation rates that are mode dependent. 

The nonstationary closed reduced equations ( |6.8| ) constitute our final result of 
the best-fit closure analysis. We emphasize that this coarse-graining of the TBH 
dynamics depends only on a single closure parameter, 7. By adjusting / the best-fit 
theory endeavors to model the linear dissipation, the modifications of the nonlinear 
interactions, and the plateau effect, all of which are mode dependent. The quantita- 
tive validation of this model via numerical experiments is described in next section. 
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7 Numerical model and experimental setup 

7.1 Direct numerical simulations (DNS) 

To ensure that the equilibrium distribution of the TBH model is close to the 



Gibbs distribution ( |2.7[ ), we choose n = 50 complex modes. This particular value 
has been tested extensively numerically in the literature ll27l . To obtain a reduced 
model that is much smaller but still resonably complicated, we choose m = 5, re- 
taining the 5 complex modes of smallest wavenumber, k = 1,...,5. Thus, our 
reduction is from 100 independent real modes to 10 real modes. 

For the prognostic numerical integration of the Fourier modes zu{t), k = 1 , . . . , 50, 
we use a fourth order Runge-Kutta timestepping method. The nonlinear advection 
terms are calculated by evaluating the Fast Fourier transform of P n (hu%\ and then 
multiplying by ik; that is, a pseudo- spectral method is utilized to integrate equa- 
tions ( |2.2| ) and (2.3 ). This algorithm has been found by ETll to respect the conser- 



vation of E to a high accuracy; throughout the time integrations reported here we 
observe conservation to an accuracy of around 10~ 5 , the expected value of E being 
10.0. 

Ensembles are constructed by sampling initial conditions from a particular trial 
distribution of the form ( |2.8[ ). The 100 real modes under these distributions are 
Gaussian with variances all equal to 1 /2j8. A value of /3 = 5.0 is adopted for the 
inverse temperature, corresponding to a standard deviation of 1 / vlO for the real 
and imaginary parts of all z.k- The means of the resolved complex modes, z\, ■ ■ . ,Z5, 
for the initial trial distribution trial distributions are derived as follows: From a 
long and hence equilibrated integration of the numerical model we draw a vector 
bk = Zk{T), k = 1, ... ,5; the initial means, ak(0), are then specified by multiply- 
ing bk by a fixed factor r</ ev . The ensemble members for the initial distribution 
are drawn by sampling the resulting Gaussian distribution. This approach allows 
us to test the sensitivity of results on the magnitude of the initial deviation from 
equilibrium. We examine the two cases, r^ ev = 1 /\/l0 and r^ ev = VlO, and refer 
to them below as the "close to" and "significantly removed from" equilibrium ex- 
periments, respectively. We comment below on other choices for r^ CT . The results 
reported here are found to be qualitatively similar when different choices for the 
bk are adopted. We thus report results from only one particular set of randomly 
generated b^. 

Since the focus of this study is on the time evolution of the resolved mean vari- 
ables, ak{t) = (z,k\p(h(t))), we construct ensembles of sufficient size so that this 
first moment is statistically steady. Emprically this has been determined to be more 
than met by samples of size 10 6 , which we adopt here. As an objective test of the 
ensemble size issue, we have examined subsample results with size 10 5 and have 
noted only very small differences to the results. As we shall discuss further below, 
the ensemble means ak(t) — > for sufficiently large t, and in fact the time scale 
for this first moment equilibration varies inversely with the wavenumber k. The 
results we report are integrated until t = 1.5, which ensure all the complex modes 
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except the first have very clearly equilibrated first moments; the first complex mode 
(k = 1) has nearly but not completely equilibrated after this integration time. 

7.2 Closure model computations 

The appropriate closure equation for comparison with the DNS results is given 



by equations (6.8 1 and (6.9 1 in Section 6. There are two important features of 



these two sets of equations. The reduced equation (6.8) controlling the first mo- 
ment evolution is similar in form to the original DNS model but with modified 
advection coefficients and with wavenumber-dependent damping. To solve it we 
modify the original TBH integrator by truncating to m complex modes and revert 
the calculation of the advection term to a purely spectral calculation that uses the 
time-dependent factors £2. These factors, which modify the advection coefficients, 
are determined by the decoupled, forced and damped, linear ODEs (|6.9[), and hence 



they are computed by a simple two time-step in which the damping and forcing are 



evaluated on the backward timestep. The wavenumber-dependent damping in (6.8 1 
is evaluated at the first time step of the fourth-order Runge-Kutta scheme for ad- 
vection. 

The closure equations have one free parameter y which we determine as fol- 
lows: The total squared difference between the 2m closure modes and the first 2m 
variables of the DNS at each time step (which is .0015) is designated the error func- 
tion for the fitting execise. It is evaluated numerically for a large number of choices 



for 7 by running both equations (6. 8 1 and (6.9 1 many times, a very cheap exercise. 
It is found to be always convex with a unique minimum which is determined by a 
manual convergence technique. 

7.3 Modified closure procedure 

The time scale separation between the retained modes and the discarded modes 
is not very sharp, a common feature for turbulent systems. One might expect in- 
tutively for this situation to adversely affect the closure performance. In order to 
evaluate this possibility we extend the closure model from 2m resolved modes to 
4m resolved modes and regard the additional 2m modes as inserting a buffer be- 
tween the fast and slow parts of the system. The extended closure model has 2m 
slow modes, 2m intermediate modes, and 16m "ignored" modes with the highest 
wavenumbers and the fastest time scales, which are regarded as a "heat bath." The 
closure model is run as before but the buffer modes are initialized at mean zero, 
while the original resolved 2m modes are initialized by r^evbk as in the previous 
subsection. 

8 Numerical results 

8.1 Close to equilibrium case 

A comparison between the first moments calculated by the DNS and those cal- 
culated using the closure model with the best choice of y = 64.74 is displayed in 



R. KLEEMAN B. TURKINGTON 




(A) Complex Mode 1 



0.01 




1.5 



(B) Complex Mode 2 

Figure 8.1. Comparison of the evolution of first moments for the di- 
rect numerical simulation (DNS) and the closure model. The initial con- 
ditions here are close to equilibrium (see text). There are five panels 
corresponding to the five retained complex modes of lowest wavenum- 
bers. In each panel the real and imaginary part of each complex mode is 
plotted for both the DNS and closure. As usual both components of the 
complex modes have the same wave number. 
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(D) Complex Mode 4 
Figure 8.1. (continued) 
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(E) Complex Mode 5 
Figure 8.1. (continued) 



Figure 8.1 Overall the RMS error per mode per timestep is found to be 0.002035, 
which may compared with the amplitude of the first moments. Each panel shows 
the behavior of the real and imaginary parts of a^{t) and all m = 5 closure modes 
and DNS modes are plotted. The performance is quite good, and in particular two 
distinctive features of the statistical dynamics are well simulated qualitatively by 
the closure: 

(1) The DNS equilibration time for each complex mode is very clearly propor- 
tional to wavenumber and excellently reproduced by the closure equations. 
This provides strong evidence that the correct dissipation for reduced mod- 
els of TBH is given by a fractional diffusion process. It is notable that 
exactly this kind of dissipation is the bare minimum required by the infi- 
nite Fourier mode Burgers-Hopf equation to prevent the appearance of a 
singular shock in finite time ifTTl . The optimal value of / reported above 
corresponds with an exponential damping time of 0.2119 /\k\, according to 
the simplified analysis at the beginning of Section 5, and this timescale is 
clearly visible in the approach to equilibrium of all modes. 

(2) Each mode exhibits essentially two well-defined regimes during equilibra- 
tion: An initial "plateau" period and a later exponential decay to equi- 
librium. These periods are of comparable duration and are both directly 
proportional to wavenumber. Both of these features are predicted by the 
theoretical analysis of Section 6 and are readily apparent in Figure |8.1| 
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most obviously for the low wavenumber cases. As the theoretical analy- 
sis shows, the fractional diffusion is ramping up during the plateau period, 
suggesting that the evolution of the resolved variables is more under the 
control of the initial conditions and to a lesser extent the non-linear inter- 
actions. 

The discrepancy between DNS and closure tends, in general, to increase with time 
and is also more apparent in the low wavenumber modes. The DNS exhibits some 
oscillatory behaviour as it approaches equilibrium, which the closure model does 
not reproduce. These oscillations have a wavelength inversely proportional to the 
mode wavenumber. Nonetheless, a careful inspection reveals that the closure tra- 
jectories tend to "bisect" these DNS oscillations cleanly in every instance. 

8.2 Significantly removed from equilibrium case 

Historically, closures theories have tended to be more successful in situations 
that are near statistical equilibrium. In the previous subsection our comparisons 
used initial mean deviations smaller than a typical equilibrium excursion of the sys- 
tem, namely, 1 / \/2p for each real mode. We now, therefore, pose a more stringent 
test of the closure to examine the case significantly removed equilibrium behavior. 
For simple expediency we multiply all deviations in the initial mean resolved vari- 
ables by a factor of \/T0, the r^ ev parameter. The initial condition amplitudes are 
therefore an order of magnitude larger than the first experiment. This amplification 
also emphasizes the contribution of the quadratic advection terms in the closure 



equations. The results are displayed in Figure 8.2 in the same format as Figure 8.1 
The best choice tunable parameter in this case is found to be y = 73.83, which cor- 
responds to an exponential damping time of 0.2602/|&|. The RMS error per mode 
and time step is 0.02561 in this instance. Since all curves have been scaled up by 
a factor of 10, this magnitude of error indicates that the fit is approximately 25% 
worse than in the close to equilibrium case. 

The general behavior of the DNS equilibration is fairly similar to the close to 
equilibrium case, but with some notable differences. In particular, the general equi- 
libration time scale is very much the same as the close to equilibrium case, again 
reinforcing the success of the reduced model with fractional diffusion. The pri- 
mary difference is in the plateau period, during which larger differences between 
the closure and the DNS are evident, particularly for the lowest two wavenum- 



bers (Figures 8.2a and 8.2b I. This deviation suggests that the modification of the 
advection terms in the closure model is somewhat less successful than the robust 
prediction of effective, fractional diffusion. 

8.3 Other initial condition deviations from equilibrium 

We also examined the case r^ ev = 1 which lies midway between the choices 
considered above. Results are qualitiatively very similar (and indeed slightly im- 
proved) over the close to equilibrium case. The best choice for damping time here 
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(B) Complex Mode 2 



Figure 8.2. Same as Figure 8.1 except the initial conditions here are 
significantly removed from equilibrium (see text) 



NONEQUILIBRIUM BURGERS DYNAMICS 



29 




(D) Complex Mode 4 
Figure 8.2. (continued) 
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(e) Complex Mode 5 
Figure 8.2. (continued) 

was 0.28156/|fc| which is little removed from the close to equlibrium case consid- 
ered above. 

In contrast the case r^ ev = 10 produces poor results. This, of course, corre- 
sponds with a very large, ten standard deviation, excursion from equilibrium. The 
RMS error per mode and time step is 0.235 and this is around four times what a 
simple rescaling of the results above would suggest. In addition the best fit damp- 
ing time was dramatically shortened to 0. 13805/|&|, about half what the previous 
three settings give. Qualitatively the results (not shown) had big discrepancies 
which were most apparent as the various modes approached equilibrium, when 
large oscillations in the DNS were not reproduced by the closure. Failure of the 
closure in this case is perhaps not a complete surprise given that a Taylor series 
expansion about equilibrium was used to solve the Hamilton- Jacobi equation (see 
Section 6). 

8.4 Robustness of tuning parameter 

The results above show that for the three settings, r^ ev = 1 / ^/l0, r^ev = 1 and 
fdev = vl0> the damping time only varies by approximately 8%. Given that the 
amplitude had been increased by an order of magnitude, the small change in the 
single adjustable parameter 7 is a rather reassuring indication of the robustness 
of the closure theory. To illustrate this robustness further we set the adjustable 
parameter 7 for the significantly removed from equilibrium case to be equal to the 
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best-fit parameter for the close to equilibrium case. The ensuing closure dynamics 
is exhibited in Figure 8.3 which should be compared with Figure 8.2 
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(D) Complex Mode 4 
Figure 8.3. (continued) 
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(E) Complex Mode 5 
Figure 8.3. (continued) 



8.5 A time scale separation experiment 

As noted previously the system under study does not exhibit a clear separation 
of time scales between the slow, low wavenumber modes and the "neglected" fast, 
high wavenumber modes. Indeed, as shown in ETl . this time scale is empirically 
observed to be inversely proportional to wave number, just like the dissipation time 
scale discussed above. Thus, for example, the second set of 2m modes have a time 
scale on average of only one third less than that of the last retained mode. This issue 
makes the construction of closure models for this system a challenge, as has been 
noted in detail in the stochastic mode reduction work of Majda, Vanden-Eijnden 
and collaborators (see Il23l ). Given this background it is reasonable to suspect that 
the discrepancies observed above, most notably for the significantly removed from 
equilibrium case, are due in part to time-scale separation issues. To test this idea 
in the context of our approach, we extend the closure model by incorporating an 
additional 2m "buffer" modes in the resolved variables, to separate the time scale 
of original closure modes from that of the neglected "heat bath" modes. In this 
new configuration the first unresolved mode has a time scale of one half of the last 
original resolved mode. 

The buffered results for significantly removed from equilibrium initial condi- 



tions may be seen in Figure 8.4 in the same format as the previous two experiments. 
Interestingly, the best fit closure now exhibits a longer damping time of 0.3004/ \k\. 
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This slower equilibration is perhaps explained by the inclusion of the buffer modes 
which mediate the nonlinear dissipation via the advection terms of the closure. 
The fit of the original modes is significantly improved relative to the unbuffered 
run, with the RMS error reduced to 0.0218. Qualitatively, this improvement occurs 
most strikingly in the higher wavenumber modes 4 and 5 , where the fit is now really 
quite impressive. This agreement might be expected since presence of the buffer 
modes removes the abrupt discontinuity in the original model at wavenumber 5. 
There are also some smaller improvements with the low wave-number modes, es- 
pecially in the plateau phase, but still the discrepancy with respect to the DNS 
oscillations on approach to equilibrium remains. 

It is interesting to compare the results above to those in 11241 where a stochas- 
tic reduction of the TBH system is proposed. These authors consider the case in 
which m = 1 , 2 and assess performance using auto-correlation rather than first mo- 
ment as is done here. In the near equilibrium case these are comparable due to the 
fluctuation dissipation relation. Like here they find that relaxation to equilibrium 
is exponential with a decay time given by the eddy turnover time. They also find 
the DNS oscillatory behaviour we observe and are also unable to produce it with 
their reduced model. They attribute this discrepancy to the lack of time separation 
between fast and slow modes. Our results with a significant buffer do not pro- 
duce a large improvement in performance of the low wavenumber variables. This 
difference deserves closer future attention. 

9 Summary and Conclusions 

A fundamental problem in nonequilibrium statistical mechanics, studies of tur- 
bulence and stochastic modelling of complex systems is the formulation of tractable 
reduced models of the slow, coherent part of the dynamical system. One of the au- 
thors has proposed a general variational principle for producing such closures in 
the context of Hamiltonian dynamics. This closure theory fits a continuous time 
series of trial probability densities to the Liouville equation using a weighted, time- 
integrated, mean-squared cost functional to quantify the lack-of-fit. The closed re- 
duced equations satisfied by minimizers of the cost functional are then obtained 
by classical Hamilton- Jacobi theory. More details may be found in the companion 
paper iPTI 

In this contribution we have applied this closure theory to a Galerkin-truncated 
Burgers-Hopf model, which has been proposed as a simple analog of more complex 
and realistic fluid systems. The model is governed by a Hamiltonian dynamics, has 
a Gibbs invariant measure, and exhibits mode decorrelation with time scales that 
vary inversely with wavenumber. These properties make it an ideal testbed for 
future work with more realistic turbulent systems. In the present work we have 
concentrated on the free relaxation of the system from nonequilibrium statistical 
initial conditions. 
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Figure 8.4. Same as Figure 8.2 except the closure model has an addi- 
tional 5 complex "buffer" modes (see text) 
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(D) Complex Mode 4 
Figure 8.4. (continued) 
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Figure 8.4. (continued) 

Our results pertain to a reduced model for the lowest 5 complex modes in a dy- 
namical system of 50 nonlinearly interacting complex modes. In light of the sim- 
plicity of the Gibbs measure, under which all the modes are independent Gaussians, 
it is possible to carry out the calculation of the governing equations for the reduced 
model explicitly up to second-order (quadratic nonlinearity) in the mean resolved 
variables. The result is a set of governing equations for the temporal evolution of 
the first moments of the resolved modes which resembles a severe truncation of 
original nonlinear system with the following modifications and characteristics. 

(1) There is a fractional dissipation which is proportional to the wavenumber 
of the mode. Remarkably this fractional diffusion is the minimum such 
regularization required to prevent finite-time singularities in the continuous 
(infinite mode) Burgers-Hopf system iTTTll . It is the primary mechanism 
by which the reduced model irreversibly equilibrates, and as expected by 
linear-response theory it has the same wavenumber-dependent time scale 
as the modal decorrelation time. 

(2) The nonlinear terms of the best-fit reduced model have modified coeffi- 
cients which imply an altered energy flux between resolved modes com- 
pared with the original TBH system. By means of these coefficients the re- 
duced model approximates the indirect influence of the unresolved modes 
on the interactions between the mean resolved modes. 
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(3) Both these modifications of the dynamics are time-dependent in that they 
take a time comparable to the equilibration time to manifest themselves in 
full. This initial period, which we term the plateau time, is incorporated 
naturally in the nonstationary version of our closed reduced model. 

The first-moment best-fit closure model has been compared with a direct nu- 
merical simulation (DNS) of the full TBH system using very large ensembles. The 
comparisons are generally good for a broad range of initial deviations from equi- 
librium. If however the deviation is made sufficiently large the fit does break down 
significantly. This is not surprising given that the closure equations rely on a Taylor 
expansion of the Hamilton- Jacobi equation centered on equilibrium. 

The wavenumber-dependent equilibration time scale predicted by the best-fit 
closure matches that produced by the DNS across the spectrum of resolved modes. 
Only one scalar parameter is adjusted empirically to fit the closure predictions 
to the DNS results. The value of this parameter is observed to be quite robust 
even as the imposed nonequilibrium initial conditions are increased by an order of 
magnitude. 

The TBH dynamics is a rather severe test for reduced models given that there 
is no clear separation of timescales. We have found that the performance of the re- 
duced model can be improved by inserting some "buffer" modes, which widens the 
separation of time scales between the slow and fast parts of the system. Nonethe- 
less, some slow oscillatory behavior of the low resolved modes is missed by the 
reduced model, indicating perhaps that there are memory effects inherent in the 
statistical dynamics. It is conceivable that an extended set of resolved variables 
that includes the time derivatives of the lowest modes might be able to capture 
those oscillations. 

Given that the TBH dynamics has been developed as a paradigm for more com- 
plex dynamical systems having quadratic nonlinearities typical of hydrodynamical 
equations of motion, the question arises whether the methodology that we have 
applied in the present investigation could be extended to other systems. When 
coarse-graining complex systems of this type onto their low, slow modes, one may 
expect to encounter Liouville residuals having similar structure to that in the TBH 
case. Under such circumstances cost functions akin to the one that we have con- 
structed in the TBH case, which quantify the collective effect of products of pairs of 
modes, may prove to be an efficacious way to represent the influence of unresolved 
fluctuations on the reduced dynamics. Future investigations on other systems of 
this kind exhibiting turbulent dynamics would be helpful in evaluating the range of 
applicability of our best-fit approach and the useful forms of the cost functions that 
define it. The reader is referred to 1131! for some further discussion on these issues. 
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Appendix: Ancilliary Material 

Here we collect some standard conventions and calculations concerning func- 
tions of several complex variables which may be useful to the reader a various 
points in the paper. 

For any smooth, complex- valued function f(z), of n complex variables, z = 
(zi, . . . ,z n ) £ C", with Zk=Xk + the usual derivatives are defined by 

d l = l ( d L_ i d L) d L = l ( d L +i d L 

dzk 2 V dx k dy k J ' dz\ 2\ dx k dy k 

the notation z\= x k — iy k is used for complex conjugate. In terms of these, the 
chain rule for the composite function f(z(t)), where t is a real variable, is simply 

±f( (t)) = V d£dZk,dl^4 = y d£dzk 
dt >> t n dz k dt + dz* k dt k t- n dzk dt ' 

The last equality uses the convention that z-k = z\, which occurs naturally in the 
context of a Fourier representation of a dynamics in terms of complex amplitudes 

Zk- 

Repeated application of this chain rule to the function f{tz), < t < 1, yields 
the complex form of Taylor's expansion, 

n n n 

f(z)=f(0)+ Uzk+ M hk 2 ZhZk 2 + N klk2k3 z kl Zk 2 Zk 3 + 0(\z\ 4 ), 

k=—n ki,k2=—n ki,k2,k3=—n 

with coefficients 

L t =g(0), »ft*-i^(0), «,m,-|^^(0). 

In Section 5 this expansion is used to determined an approximation to the value 
function, v(A), for A G C m , and in Section 6 for the time-dependent value func- 
tion v(X,t). For those calculations this form of the expansion is better than the 
equivalent multi-index form. 

The Lagrangian if (A, A) and the Hamiltonian J^(A,/i) that occur in the defin- 
ing optimization principle and the associated Hamilton- Jacobi equation, respec- 
tively, are related by the Legendre tranform (4.8 1. It is instructive to relate this 
complexified Legendre transform to the familiar real transform. Write the com- 
plex variables as X k = ^ k + iT] k and \i k = <p k + i\j/h for & = 1 , . . . , m. Then, (4.8 1 is 
equivalent to 

d if 5 if ™ , JSf 

d£, k 2 dT] k 2 2 2 



Thus the complex Legendre transform ( |4.8| > between if and is identical with 
the r eal tranform between if/2 and Jt/2. The stationary value function, v(A), in 
(3.8 1 is defined by the action integral for if (A, A) and satisfies the complex form 
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of the Hamilton- Jacobi equation (4.9 1 
by v is 



The conjugacy relation (|4.10[) determined 



_d_ ._d_ 
d& dr\ k 



This relation shows that ( |4.9| > and ( |4.10[ ) are identical with the real Hamilton- 
Jacobi theory applied to the Hamiltonian , T] , , / 2 and the value function 
v(^,tj)/2. The nonstationary case is complexified in the same way. The common 
factor of 1/2 throughout is irrelevant to the analysis, being absorbed in the complex 
variable conventions. 
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